Computational modelling of the collective stochastic motion of Kinesin nano motors 
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We have developed a two dimensional stochastic molecular dynamics model for the description of 
intra cellular collective motion of bio motors, in particular Kinesins, on a microtubular track. The 
model is capable or reproducing the hand-over-hand mechanism of the directed motion along the 
microtubule. The model gives the average directed velocity and the current of Kinesins along the 
microtubule. It is shown that beyond a certain density of Kinesins, the average velocity and current 
undergo notable decrease which is due to formation of traffic jams in the system. 



I. INTRODUCTION 



Various complex functions of Eucariotic cells namely 
mitosis, intracelluler transports and cell motility are 
based on the cooperative motion of molecular motorpro- 
teins such as Dyneins, Kinesins and Myosins P, 0,11,11] ■ 
These molecular motors move on cytoskeletal filaments 
like microtubles and F-actins. The mentioned phenom- 
ena exhibit a variety of self-organized processes and 
patterns which are characterized by their time scales 
E i, 0, B A large number of in-vitro experi- 

ments have revealed that the motion of Kinesin is pro- 
cessive i.e., it performs stroke- type steps on its track 
before detachment [l^, [HI, [13 ■ It is now well estab- 
lished that the Kinesins motion is derived by the free en- 
ergy released in the chemical hydrolysis reaction of ATP 
(adenosine-triphosphate). The fuel of the motion that 
is the ATP molecules arrive randomly in the vicinity of 
Kinesins. Consequently, the directed motion takes place 
on a stochastic grounds. During each step in which one 
period of the microtubule, hereafter referred to as MT, 
is covered, one ATP molecule is hydrolysed. The under- 
lying mechanism responsible for this forward motion is 
partially understood. Two basic mechanisms have been 
proposed. In the first one the so-called inchworm mech- 
anism, one head of Kinesin drag the other one which 
is always at the rear [Tsj . In this mechanism, there is 
an asymmetry between the lagging. The second mech- 
anism which is called hand-over-hand has much resem- 
blance to human walking. This mechanism involves re- 
peated docking and undocking of both heads upon break- 
ing of the chemical bond with the MT in the course of 
ATP hydrolysis At each stepping event, the rear 

head executes a forward motion of 16 nm and attaches 
to the tubulin binding site ahead of the front head. Re- 
cent empirical evidences overwhelmingly have established 
that hand- over- hand mechanism is the true mechanism 
p^ . In recent years various approaches have been pro- 
posed to model the directed motion of nano-sized bio mo- 
tors. Simple random walk was used as a simple though 
efficient framework for modelling such kind of motion 

[as 

, [TSi . ,19,] . In this framework, Kinesins are assumed 
to be point particles which execute a directed random 
walk on a discrete one dimensional filament, as a model 



for microtubular track, with lattice spacing equal to Ki- 
nesin walking distance i.e. 8 nm. Due to stochasticity in 
the Kinesin motion, it would be reasonable to incorpo- 
rate degrees of randomness in the modelling approach. In 
this spirit, a new approach namely asymmetric particle 
hopping models in continuous time which are based on 
master equation approach came into play [13, [13, [HI [3 ■ 
Another theoretical approach for the modelling of bio mo- 
tors motion is the so- ca/ted ratchet mechanism R|. In this 
mechanism, the directed motion is generated by a time- 
dependent periodic potential which is switched between 
two states [m [H, [11, [H, [H [H, [H [s^l • The stochastic- 
ity comes in the manner in which the potential switches 
from one state to the other one and reflects the random 
arrival of ATP molecules and conformational transitions 
which changes the interacting potential between Kinesin 
and MT. In principle, more than one Kinesins can simul- 
taneously move on a single MT. Kinesins do certainly 
interact with each other. While the true interaction po- 
tential between adjacent Kinesins has not yet empirically 
explored, in-vitro observations have revealed the collec- 
tive motion of Kinesins which arises from the interaction 
between these motors fsl]. In order to find some insight 
into the nature of the Kinesins motion, we have devel- 
oped a Langevin-type model which incorporates interac- 
tion among Kinesins. Our aim is to improve our under- 
standing on the spatio-temporal organisations of intra 
cellular transport phenomena and find the mechanism 
responsible for formation of molecular traffic jams. 



II. MODELLING OF KINESINS MOVEMENT 

The model we have considered is a generalisation of the 
model proposed in |30j for the motion of a single Kinesin. 

We shall now explain our stochastic ratchet model in 
some details. The force acting on each Kinesin head con- 
sists of four terms as follows: 
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Hratchet IS the stochastic ratchet potential which mod- 
els the interaction between Kinesin and microtubule. 
Hbistabie IS a bistable potential that models the elastic 
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coupling between two heads of Kinesin, H^ep is the re- 
pulsive potential between adjacent Kinesins and finally 
r represents both the stochastic Brownian and the fric- 
tional forces experienced by Kinesin heads due to thermal 
noise in the intra cellular environment. Let us explain 
these four terms in more details. 



A. Ratchet potential 

The ratchet potential H{x, y, t) should involve an 
asymmetry along the MT. More precisely, we have taken 
the form of ratchet potential as follows [ 30l |: 



H(x,y,t) = CH,{x)Hy{y)A{t)S{t) 



(2) 



With periodicity in x direction i.e., H{x + L,y,t) — 
H[x, y, t) where L is the microtubule period which is 
about 8 nm. For Hx{x) we have used the first eleven 
terms of a Fourier series of an asymmetric function W{x) 
as follows: 



.r-^ 2'Kmx. , . 2Trm.x 

Hx(x) ^ ao+ arnCos{ — - — ) + hmSin{ ^ ) (3) 



The coefficients ao,am and hm m — I,-- - ,5 are 
Fourier coefficients. For the asymmetric function W{x) 
we chose the saw-tooth function: 



w^(x) = io(| - [|]) - 10 I 



> 0.9 



(4) 



Where [ ] denotes the integer part. The summation 
form of Hx{x) has been used for computational conve- 
nience. Figure (1) shows Hx{x) versus x . 

Each peak of Hx{x) is located at a distance O.IL from 
the position of the next minimum. For more details 
see Ref. [131 . Concerning the part of potential Hy{y) 
responsible for the perpendicular direction to the MT 
we have used the following potential : 



Hy{y) = - 



y- 



(5) 



where yo = 0-1 Vi. It is assumed that the above form 
of Hy (y) can model the interaction between the molecules 
making up the Kinesin head and the MT. The symme- 
try of the function allows for detaching both in upward 
and downward directions perpendicular to the MT. This 
function is smoother than the Morse potential therefore 
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FIG. 1: Periodic dependence of Hx on x. 




FIG. 2: Dependence of Hy on the perpendicular direction y. 



we can choose larger time steps which can speed up the 
simulation. Figure (2) shows Hy. 

Now we shall discuss the time dependent terms i.e., 
A{t) and S{t) of the Ratchet potential. The stochas- 
tic variable A{t) specifies the biochemical affinity of each 
Kinesin head. It is a flashing potential which can switch 
between two discrete values A{t) = 0, 1. The value of 
A{t) depends on the biochemical status of the ATPase 
binding site of the head. The biochemical cycle is divided 
into four status. [13, [13 In the first status, denoted by 
K, the ATPase binding site contains no ATP molecule. 
When an ATP molecules arrives at the binding site, the 
head is said to be in its second status which is denoted by 
K.ATP. The third status is associated with the hydrolysis 
of the absorbed ATP into ADP plus an inorganic phos- 
phate molecule Pi. This status is denoted by K.ADP.Pi. 
Eventually after releasing the Pi molecule, the head is 
characterized by its fourth status K.ADP. The biochemi- 
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cal cycle is accomplished by releasing the ADP molecule. 
In the status 1-3 the head has a high affinity {A = 1) 
to attach, or to keep its attached condition to the MT. 
On the other hand, in the fourth status i.e., K.ADP the 
affinity will considerably reduced {A = 0) due to some 
conformational changes. The second term S{t) is related 
to the interaction between Kinesins. Besides biochemical 
status, each head can be either attached to the MT or 
detached from it. Each head can attach only to partic- 
ular sites on the MT. Since only one head can attach to 
each MT binding site, the potential seen by each head 
depends on the occupancy of its adjacent MT binding 
site. For adjacent Kinesins to MT, the stochastic vari- 
able S{t) describes the occupancy of the corresponding 
nearest binding site on the MT. It is zero if this site is 
already occupied by another head and one otherwise (see 
Fig. 4 ). Parameter C is set at C = 0.01 eV. The value of 
constant C determines the depth of the potential well. Its 
size must be such that the energy difference between the 
on-state (A=l) and the off-state (A=0) under the equi- 
librium condition is comparable with the energy released 
from the hydrolysis of the ATP i.e., 

[H{x,y,t)]A{t)=i - [H{x,y,t)]A{t)=o < EATP^ADP+Pi 

(6) 

B. Viscous and stochastic thermal noise terms 

The force F, which for simplicity acts in two dimension, 
has the following structure: 

=-?7i(t)+ex(i). (7) 

Ty^-j^y{t)+iy{t). (8) 

The first term corresponds to energy dissipation due 
to viscosity of the surrounding fluid. and £,y are sta- 
tionary random Gaussian white noises with zero mean. 
More precisely we have: 

{Ut)Ui'))='^vkBT5{t-t')5,, (9) 

The drag coefficient 77 for a Kinesin in an aqueous en- 
vironment is reported at 77 « 6 x 10~^pNs/nm We 
have taken this value for rj in this paper. 

C. The bistable potential 

The bistable potential is responsible for the elastic cou- 
pling of the two Kinesin heads. This potential is ex- 
pressed by the following form [H, HO] : 

Hh^stableiM = Cl[l + {^Y - 2(^)2] (10) 
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FIG. 3: Bistable potential between two heads vs distance. It 
models the elastic interaction between Kinesin heads. 

In which Ar is the distance between the two heads, 
Ci = 0.12 eV is the amplitude of the potential and repre- 
sents the coupling strength between heads, and I — 0.7 L 
is the distance between the two potential minima. In the 
following figure, we have sketched the dependence of the 
bistable potential versus distance. 

D. The potential between different Kinesins heads 

In order to model the interaction between adjacent Ki- 
nesins, we introduce a short-ranged repulsive potential 
between heads of different Kinesins. More concisely, we 
have used the shifted-force repulsive part of the Lennard- 
Jones potential as follows: 

HrepiAr) = f/(Ar) - U{Ar,) (11) 

In which [/(Ar) = with e = 0.2 eV, a = 70 A° 

and the cut off length Tc = 80 A". This repulsive poten- 
tial prevents the heads of adjacent Kinesins from getting 
too close to each other. 



E. simulation details 

We have not yet discussed in details how quantities A 
an S vary with time. To this end, we introduce a quantity 
a(t) for each Kinesin's head which specifies whether at 
time t the head is attached to the MT or not. a{t) is 
one if the head is attached to a binding site and zero 
otherwise. Before proceeding further let us specify when 
we consider a head attached to a MT binding site. If at 
any time step, a head has a high affinity {A{t) = 1) and 
is located within a rectangular box with size Sx = 3.8 nm 
and 6y = 0.8 nm centred at a binding site of the MT then 
we consider this head as an attached one. 
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FIG. 4: For adjacent heads to MT, The ratchet potential is 
zero if the nearest binding site on the MT is aheady occupied 
by another head. A head is considered attached to a binding 
site if it is in a smaU vicinity of that site. 



Not we turn into A{t). Apparently the biochemical 
cycles of each head can be described by transition rates 
from one state to the other one. We assume that the av- 
erage time that each head spends in the states K, K.ATP 
or K.ADP.Pi is in the order of ms when the head is at- 
tached to the MT [33, [11]. We note that this average time 
depend on the attachment of the head and the position of 
its partner head. The average time a head spends in the 
fourth state K.ADP is much shorter and is in the order 
of fis [sil, [3^ . Since the rate of ATP consumption is con- 
siderably reduced when a head is detached [13, [H, [1^ , 
in our model we have assumed that only an attached 
head received an ATP molecule. Moreover, we experi- 
mentally know the rate at which an attached head makes 
a transition from the state K — > K.ATP strongly depends 
on the attaching state and position of its partner head 
[39I [ioj . In case the partner head is attached and rear, 
the above rate is lowered about two orders of magnitude. 
We have tuned the corresponding rates in our code such 
that the processivity and the average velocity of a sin- 
gle Kinesin coincides with the known empirical findings. 
Consequently, we used these adjusted rates in the col- 
lective motion of Kinesins. Periodic boundary condition 
have been imposed and the Kinesins motion are assumed, 
for simplicity, to be restricted to two dimension along x 
and y axes. The size of the simulation box was taken as 
follows: L,, = 50 L, Ly = 40 L. The MT is located along 
the X— axis at y = and y £ [— j^, -j-]- Time step St 
was set to 5 X 10~^°s. The system has been simulated for 
3 X 10^° timesteps in each run. The first 5 x 10^ timesteps 
are discarded for reaching to equilibrium. Kinesin global 
density pg is defined as the number of Kinesins per bind- 
ing site of the MT. Certainly pg can exceed one. 



III. SIMULATION RESULTS 

Generically the Kinesin movement can be classified 
into 3 states: resting, forward proccesive movement and 
detachent/atachement to the MT. Once a Kinesin leaves 
the box from the vertical boundaries at y = its cor- 
responding image re enters from the opposite side. Those 
Kinesins which leave the last site in the x— direction at 
i = Lx enter into the MT from the left site i = 1. 
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FIG. 5: 2D trajectory of a typical Kinesin at pg — 0.4. Its 
corresponding space-time plot is sketched in Fig. 7. 
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FIG. 6: Space-time plots of some Kinesins and the time de- 
pendence of their centre of mass at pg = 0.4. Despite each Ki- 
nesin undergoes rapid fluctuations upon becoming unbound, 
the centre of mass (CM) has a smooth increasing behaviour. 



Figure (5) exhibits the trajectory of a typical Kinesin at 
the global density pg = 0.4: 

Figure (6) exhibits the space-time plots of some more 
Kinesins and the time dependence of their centre of mass 
at Pg = 0.4. We observe that although each Kinesin un- 
dergoes rapid fluctuations upon becoming unbound, the 
centre of mass (CM) has a smooth increasing behaviour 
in time. We have found the CM velocity by fitting a linear 
line to the the curve as Vcm = 0.5 pm/s. As can be seen, 
when a Kinesin is attached to the MT, its x component 
undergoes small fluctuations and is frequently increasing 
with time. It rarely undergoes backward motion which 
is consistent with in- vitro experiments [1, [4l[, \^ . 

On the occasion of detachment from the MT, the Ki- 
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FIG. 7: Processive motion of a Kinesin along the MT. 

nesin performs a stochastic random motion, in the cy- 
toplasmic environment, which is characterized by large 
fluctuations. It is a well-known fact that Kinesins are 
processive motors i.e., they can move dircctcdly along 
the filament for relatively large distances before detach- 
ing from it. Our model can successfully reproduce this 
behaviour in the collective motion of Kinesins. To illus- 
trate this behaviour, we have explicitly shown x compo- 
nent time dependence of the two heads as well as the 
centre of mass for a typical Kinesin in the following fig- 
ure. 

As can be seen from Figure (7), the hand-over-hand 
mechanism of the directed motion is evident. In order to 
find a better insight, we have evaluated the dependence 
of a single Kinesisn mean squared displacement (MSD) 
on time. The MSD has been obtained by averaging over 
the trajectories of all the Kinesins. Figure (8) depicts 
the time behaviour of a single Kinesin' MSD for various 
global Kinesin concentration pg in a log-log plot. 

We have evaluated the slopes by linear curve fitting. 
To a very good approximation the diffusion is normal 
for densities larger than pg ~ 0.2 i.e., the MSD grows 
linearly in time. It would be more appropriate to define 
a line density p. This quantity is defined as the time 
average number of bound Kinesins to the MT. Note that 
a bound Kinesin can occupy two or one site of the MT 
depending on whether one or two heads are attached. We 
have sketched the dependence of p on pg in Figure (9). 

We observe a linear increase of p up to pg = 0.4. After- 
wards the rate of increase is lowered. This is due strong 
repulsion between bound Kinesins. We stress that in or- 
der to have values of p larger than 0.6 we have to sig- 
nificantly increase Pg. This enormously rises the com- 
putational costs. We extrapolated the behaviour of the 
curve for larger p. It turns out that p tends to an asymp- 
totic value near 0.7. More specifically, let us denote the 
number of those bound Kinesins attached via one head 
to MT by Ni and those attached via two heads by N2 so 



FIG. 8: Single Kinesin mean squared displacement versus 
time for various densities. To a very good approximation the 
diffusion is normal for densities larger than pg ~ 0.2 i.e.; the 
MSD grows linearly in time. 
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FIG. 9: Dependence of the line density p on the global density 
Pg. p tends to an asymptotic value near 0.7. 

we have N = Ni) + Ni + N2 (N is the total number of 
Kinesins). Furthermore, we show the number unbound 
Kinesins by Nq and the MT binding sites by Nbs- 

N0 + N1+N2 N1 + N2 

= Ws ' ^ = ^v^ ^''^ 

In the extreme limit Pg — ^ 00 we have (by extrapola- 
tion) p = 0.7, on the other hand we have 2N2 + N1 = Nbs 
which means all the sites arc occupied by Kinesin heads. 
We therefore find that sixty percent of the sites are oc- 
cupied by double-head attached Kinesins and forty per- 
cent by single-head attached Kinesins. In the following 
graphs, we exhibit the dependence of some cooperative 
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FIG. 10: Diffusion coefficient D of a single Kinesin versus den- 
sity. The diffusion constant D turns out to be an increasing 
function of p. It grows slowly for small densities but shows a 
rapid increment for p > 0.3. The reason is due to increase of 
unbound Kinesins when the density is enhanced. 



FIG. 11: Averaged CM velocity versus density p. (T4) is 
a decreasing function of p. For small densities it decreases 
smoothly. When p goes beyond p = 0.2, {14) shows a rapid 
decrease. 



quantities on the line density rho. In Figure (10) we 
exhibit the dependence of a single Kinesin diffusion co- 
efficient D versus p. It grows slowly for small densities 
but shows a rapid increment for p > 0.3. The reason is 
due to increase of unbound Kinesins when the density 
is enhanced. Obviously, an unbound Kinesin performs a 
diffusive motion therefore increasing the number of un- 
bound Kinesins gives rise to enhancement of the overall 
diffusion which in turn leads to a larger D. Note that D 
has been obtained by averaging over the trajectories of 
all the Kinesins. 

We next turn into the issue of directed motion of 
Kinesins. First, we have computed the dependence of 
the averaged velocity of the CM versus p. Figure (11) 
sketches this behaviour. 

To obtain the above graph, we plotted Xcm versus time 
and fitted a linear curve to it. The slope of the fitted 
line corresponds to the averaged velocity. As expected, 
{Vx) is a decreasing function of p. For small densities 
it decreases smoothly. However, when the density goes 
beyond p = 0.2, {Vx) shows a rapid decrease. In figure 
(12) we have exhibited the mean directed passage length 
(processive run length) , the mean processive time and the 
velocity of this directed motion on the MT as a function 
of p. 

The average processive run length Lprs defined as the 
average distance a Kinesin can move proccesively on the 
MT before detachment exhibits a rapidly decreasing be- 
haviour up to p = 0.4. After this value it shows a weak 
dependence on p. We speculate that at this density large 
traffic jams are formed which do not allow Kinesins to 
move directedly. The analogous behaviour is observed 
for the mean temporal processivity Tp^s- Dividing Lp^-s 
by Tprs gives us the average velocity Vprs of the proces- 
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FIG. 12: Averaged directed processive run length, time and 
velocity versus density. Average processive run length ex- 
hibits a rapidly decreasing behaviour up to p = 0.4. After 
this value it shows a weak dependence on p. The analogous 
behaviour is observed for the mean temporal processivity. 



sive motion. The average velocity has a more smooth de- 
creasing behaviour. Having obtained the velocity of the 
directed motion, we are able to find the dependence of the 
current of Kinesins along the microtubule. The current 
J is defined by the fluid mechanics relation J — PgVprs- 
Figure (13) sketches the dependence of J versus p. 

J exhibits a maximum at p ~ 0.34. The interesting 
point is the existence of an asymmetry in the J — p di- 
agram. Normally in lattice driven gas models such as 
asymmetric simple exclusion process (ASEP) J appears 
as a symmetric function of p [43l |. The current J is re- 
lated to the rate of cargo transport by Kinesin motor pro- 
teins. In order to flnd a deeper understanding, we have 
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FIG. 13: Average Kinesins current versus density. It shows a 
niciximum at p ^ 0.34. 



evaluated the steady state percentage of the bound and 
unbound Kinesins as a function of p. Bound Kinesins are 
divided into two groups: one-head attached and two-head 
attached to the MT. Figure (14) exhibits this behaviour. 
Percentage of unbound Kinesins increases with density. 
In the same manner, percentage of one-head attached Ki- 
nesins increases with p but it seems to become saturated 
at large densities. Interestingly, percentage of two-head 
attached Kinesins shows a decreasing behaviour. The 
rate of decrease is small at low densities but becomes 
larger for higher densities. 

IV. SUMMARY AND CONCLUDING 
REMARKS 



100 



80 



20 







1 1 1 1 1 1 1 




A 


- Detached kinesin \ 




— a — 


- 1 head atached kinesin V 




— 


- 2 heads attached kinesin - 



"o 0.2 0.4 0.6 



P 

FIG. 14: Percentage of averaged bound and unbound Kinesins 
on the MT vs density. 



We have simulated the collective motion of Kinesins 
on a microtubular track by developing a two dimensional 
Langevin-type model. The model is capable of repro- 
ducing the hand- over- hand mechanism of the directed 
motion along the microtubule. Various quantities such 
as difFiision constant, average velocity and the number 
current of Kinesins have been obtained by performing 
extensive molecular dynamics simulations. Dependence 
of the above quantities on the the overall density of Ki- 
nesins have been obtained. It is shown that beyond a 
certain density, the average velocity and the current un- 
dergo notable decrease which is due to formation of traffic 
jams in the system. This observation is in accord to the 
ubiquitous flow-density characteristics. 
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